;+
; NAME:
;       FILL_IN_RB3_MORPH
;
; PURPOSE:
;       Fill in the rb3 structure with relevant morphological matching
;       information needed for later analyses.
;
; CATEGORY:
;       IRDC Distance Analysis
;
; CALLING SEQUENCE:
;       FILL_IN_RB3_MORPH
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       CASENO -- Case number to use.  [Default: 4].
;       SUFF   -- File suffix to use (i.e. '_q8').  [Default: '']
;
; KEYWORD PARAMETERS:
;       CLEM   -- Use the Clemens (1985) rotation curve for
;                 calculation of model fit distances from N/F fit
;                 placements.
;
; OUTPUTS:
;       NONE
;
; OPTIONAL OUTPUTS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  07/04/12, TPEB -- Code from check_grs_match.pro,
;                                   consolidated into a routine with
;                                   optional inputs.
;
;-

FUNCTION DPDF_ERROR_BAR, level, TOTPROB=totprob
  
  ind = WHERE(totprob GE level[0], nind)
  area = (nind EQ 0) ? 0.d : TOTAL(totprob[ind])
  
  ;; print,'Returning: ',abs(area/TOTAL(totprob) - 0.6826d)
  
  RETURN, abs(area/TOTAL(totprob) - 0.6826d)
END

PRO FILL_IN_RB3_MORPH, caseno, CLEM=clem, SUFF=suff, DRIVEMORPH=drivemorph, $
                       ALPHA=alpha
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Parse inputs
  IF n_elements(caseno) EQ 0 THEN caseno=2
  IF caseno LT 0 || caseno GT 3 THEN message,'Improper model case number!'
  cno = string(caseno,format="(I1)")
  rotc = KEYWORD_SET(clem) ? 2b : 0b
  IF n_elements(suff) EQ 0 THEN suff = ''
  drivemorph = KEYWORD_SET(drivemorph)
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART
  defsysv, '!IRDC', exists = exists
  IF ~ exists THEN irdc_params
  
  ;; Load everything
  IF n_elements(rb3) EQ 0 THEN $
     restore,'irdc_dist_model/bgps_rb3'+suff+'.sav',/ver
  n = n_elements(rb3)
  
  ;;Check that FILL_IN_RB3_GRS has been run already...
  IF total(rb3.grs) EQ 0 THEN fill_in_rb3_grs,caseno,CLEM=clem,SUFF=suff
  
  IF drivemorph THEN BEGIN
     fmeas = 1.d - rb3.c_meas / $
             (1.d - exp(-rb3.s_peak * flux2tau_bgps(!IRDC.Td, /mjy)))
     driver_morph,cnum_list=rb3.cnum,SUFF=suff,/save,CLEM=clem,tophat=0.,$
                                ALPHA=alpha;,$
                                ; /BEAMMASK,FMEAS=fmeas,MDIST=rb3.fitdist
  ENDIF

  IF n_elements(alpha) EQ 0 THEN BEGIN
     alpha = 2.
     suff_alpha = ''
  ENDIF ELSE suff_alpha = string(alpha,format="('_',F0.2)")
  
  ;; Load remaining SAVE files
  restore,'./irdc_dist_model/irdc_morph_pvec'+suff_alpha+'.sav',/ver
  restore,'./irdc_dist_model/data/case'+cno+'_irdc_TdFIT'+suff+'.sav',/ver
  restore,'./local/bgps_velocity_struct.sav',/ver
  s = read_bgps_csv()
  
  ;; Create arrays to be used here...
  tag     = bytarr(n)
  
  ;; Go through each source and do the KDIST
  d = dindgen(!MW.NBINS)*!MW.BINSIZE + !MW.BINSTART
  FOR i=0L, n-1 DO BEGIN
     IF (i+1) MOD 100 EQ 0 THEN $
        print,'Processing object '+string(i+1,format="(I3)")+$
              ' of '+string(n,format="(I0)")
     
     this  = WHERE(velocity.cnum EQ rb3[i].cnum, nind)
     kprob = kdist_spectrum(rb3[i].l, rb3[i].b, velocity[this].v_std, $
                            velocity[this].densegas, CLEM=clem)
     
     tandist = !MW.R0*cos(rb3[i].l * !dtor)/cos(rb3[i].b * !dtor)/1.d3
     
     
     
     ;; INVX
     totprob = kprob * pvec[this].invx * prob_lonlat(s[this])
     totprob /= max(totprob,ind)
     tag[i] = (d[ind]/1.d3 LE tandist) ? 0 : 1
     rb3[i].morph_d = d[ind]/1.d3
     dind    = WHERE(d/1.d3 LE tandist, ndind)
     find    = WHERE(d/1.d3 GT tandist, nfinf)
     rb3[i].morph_np = TOTAL(totprob[dind]) / TOTAL(totprob)
     
     IF (rb3[i].morph_np GE 0.45) && (rb3[i].morph_np LE 0.55) THEN $
        rb3[i].equal_prob = 1b
     
     ;; ;; Analyze Prior Prob for width...
     ;; ;;print,size(d,/dim),size(pvec[this].invx,/dim)
     ;; yfit = MPFITPEAK(d, pvec[this].invx, A, NTERMS=4)
     ;; rb3[i].morph_width = A[2]
     
     ;; ;; Try FWHM
     ;; hm = 1.d - ((1.d - A[3]) / 2.)
     ;; hind = WHERE(pvec[this].invx GE hm, nhm)
     ;; fwind = long(minmax(hind))
     ;; rb3[i].morph_fwhm = d[fwind[1]] - d[fwind[0]]
     
     ;; Want information about the posterior PDF
     parinfo = {Value:0.5d,limited:[1b,1b],limits:[0.d,1.d],$
               tnside:0,parname:'Level',step:0.05}
     functargs = {TOTPROB:totprob}
     level = TNMIN('DPDF_ERROR_BAR',FUNCTARGS=functargs,/AUTODERIVATIVE,$
                   PARINFO=parinfo,STATUS=status,ERRMSG=errmsg,/QUIET)
     
     ;; print,'Level: ',level
     ;; plot,d,totprob,psym=2
     ;; vline,/horiz,level,color=cgColor('Yellow')
     
     ;; IF level EQ 0.5 THEN STOP
     rind = WHERE(totprob GE level, nri)
     
     print,'Area: ',total(totprob[rind]/total(totprob))
     ebars = [d[ind] - min(d[rind]), max(d[rind])-d[ind]]
     
     ;; print,'TNMIN Output: ',status,'  ',errmsg
     print,'Error Bars: ',ebars/1.d3
     
     rb3[i].morph_sigm = ebars[0]/1.d3
     rb3[i].morph_sigp = ebars[1]/1.d3
          
     ;; Ratio of peaks
     peak_np = max(totprob[dind])
     peak_fp = max(totprob[find])
     
     rb3[i].morph_pr = (peak_np GE peak_fp) ? peak_np/peak_fp : peak_fp/peak_np
     
     
     kdn = KDIST(rb3[i].l, rb3[i].b, rb3[i].vlsr, /NEAR)/1.d3
     kdf = KDIST(rb3[i].l, rb3[i].b, rb3[i].vlsr, /FAR)/1.d3
     
     ;; Determine whether an object is well-constrained,
     ;;   morphologically speaking...
     chp = (rb3[i].morph_np LT 0.5) ? 1.-rb3[i].morph_np : rb3[i].morph_np
     
     IF (rb3[i].morph_pr GE 2.) OR (chp GE 0.68) OR (abs(kdf-kdn) LE 1.0) THEN $
        rb3[i].morph_constrain = 1b
     
  ENDFOR  ;; Loop for sources
  
  print,m4_stat(tag)
  
  ;; Save rb3 file again
  save,rb3,filename='./irdc_dist_model/bgps_rb3'+suff+'.sav',/ver
  
END
